clc;
clear;

idx1 = 8;
idx2 = 4;
syms x y t nu
switch idx1
    case 1
        u1 = x.^2.*y.^2+exp(-y);
        u2 = -(2/3)*x.*y.^3+2-pi*sin(pi*x);
        p = -(2-pi*sin(pi*x)).*cos(2*pi*y);
    case 2
        u1 = (x.^2.*y.^2+exp(-y)).*cos(2*pi*t);
        u2 = (-(2/3)*x.*y.^3+2-pi*sin(pi*x)).*cos(2*pi*t);
        p = -(2-pi*sin(pi*x)).*cos(2*pi*y).*cos(2*pi*t);
    case 3
        u1 = pi*sin(2*pi*y).*sin(pi*x).^2;
        u2 = -pi*sin(2*pi*x).*sin(pi*y).^2;
        p = 10*cos(pi*x).*sin(pi*y);
    case 4
        u1 = pi*sin(2*pi*y).*sin(pi*x).^2.*sin(t);
        u2 = -pi*sin(2*pi*x).*sin(pi*y).^2.*sin(t);
        p = 10*cos(pi*x).*sin(pi*y).*sin(t);
    case 5
        u1 = 10.*x.^2.*y.*(2.*y-1).*(x-1).^2.*(y-1);
        u2 = -10.*x.*y.^2.*(2.*x-1).*(x-1).*(y-1).^2;
        p = (20.*x-10).*(2.*y-1);
    case 6
        u1 = 10.*x.^2.*y.*(2.*y-1).*(x-1).^2.*(y-1).*cos(t);
        u2 = -10.*x.*y.^2.*(2.*x-1).*(x-1).*(y-1).^2.*cos(t);
        p = (20.*x-10).*(2.*y-1).*cos(t);
    case 7
        u1 = 5*(x.^2-2*x.^3+x.^4).*(2*y-6*y.^2+4*y.^3);
        u2 = -5*(y.^2-2*y.^3+y.^4).*(2*x-6*x.^2+4*x.^3);
        p = 5*(2*x-1).*(2*y-1);
    case 8
        u1 = 5*(x.^2-2*x.^3+x.^4).*(2*y-6*y.^2+4*y.^3).*cos(2*pi*t);
        u2 = -5*(y.^2-2*y.^3+y.^4).*(2*x-6*x.^2+4*x.^3).*cos(2*pi*t);
        p = 5*(2*x-1).*(2*y-1).*cos(2*pi*t);
end
u1t = diff(u1,t);
u1x = diff(u1,x);
u1y = diff(u1,y);
u1xx = diff(u1x,x);
u1yy = diff(u1y,y);
u1dx = matlabFunction(simplify(u1x), "Vars", {x,y,t})
u1dy = matlabFunction(simplify(u1y), "Vars", {x,y,t})
u2t = diff(u2,t);
u2x = diff(u2,x);
u2y = diff(u2,y);
u2xx = diff(u2x,x);
u2yy = diff(u2y,y);
u2dx = matlabFunction(simplify(u2x), "Vars", {x,y,t})
u2dy = matlabFunction(simplify(u2y), "Vars", {x,y,t})
px = diff(p,x);
py = diff(p,y);
pdx = matlabFunction(simplify(px), "Vars", {x,y,t})
pdy = matlabFunction(simplify(py), "Vars", {x,y,t})
% u2xy = diff(u2x,y);
% u1yx = diff(u1y,x);
switch idx2
    case 1
        disp("定常stokes:");
        f1 = simplify(-nu*(u1xx + u1yy) + px);
        f2 = simplify(-nu*(u2xx + u2yy) + py);
    case 2
        disp("定常NS:");
        f1 = simplify(u1*u1x + u2*u1y - nu*(u1xx + u1yy) + px);
        f2 = simplify(u1*u2x + u2*u2y - nu*(u2xx + u2yy) + py);
    case 3
        disp("非定常stokes:");
        f1 = simplify(u1t - nu*(u1xx + u1yy) + px);
        f2 = simplify(u2t - nu*(u2xx + u2yy) + py);
    case 4
        disp("非定常NS:");
        f1 = simplify(u1t + u1*u1x + u2*u1y - nu*(u1xx + u1yy) + px);
        f2 = simplify(u2t + u1*u2x + u2*u2y - nu*(u2xx + u2yy) + py);
end
f1 = matlabFunction(f1, "Vars", {x,y,t,nu})
f2 = matlabFunction(f2, "Vars", {x,y,t,nu})